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We study the quantum phase transition between the superfluid and the Mott insulator in the 
one-dimensional (ID) Bose-Hubbard model. Using the time-evolving block decimation method, we 
numerically calculate the tunneling splitting of two macroscopically distinct states with difi'erent 
winding numbers. From the scaling of the tunneling splitting with respect to the system size, 
we determine the critical point of the superfluid to Mott insulator transition for arbitrary integer 
filling factors. We find that the critical values versus the filling factor in ID, 2D, and 3D are well 
approximated by a simple analytical function. We also discuss the condition for determining the 
transition point from a perspective of the instanton method. 
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I. INTRODUCTION 

I Systems of cold atoms in optical lattices have provided 
a highly controllable testing ground for quantum many- 
body physics Particularly, a transition from super- 
fluid (SF) to Mott-insulator (MI) has attracted much at- 
tention. The SF to MI transition can be induced by 
increasing the depth of the optical latticepotential, and 
has been experimentally realized in ID @t3, 2D [6|-Q, 

; and 3D 

It has been established that a system of cold bosonic 
I atoms in an optical lattice can be quantitatively de- 
'. scribed by the Bose-Hubbard (BH) model 

L L 

■ H=-jY.{byb,+,+h.c)+-J2n,{h,~l), (1) 

when the lattice is sufficiently deep, i.e. in the tight bind- 
ing regime. Here bj annihilates a boson at the lowest 
level localized on the jth site and fij is the number oper- 
ator. The hopping energy J corresponds to the overlap 
integral of two nearest-neighboring Wannier orbitals, and 
] decreases exponentially when the lattice depth increases. 
The onsite interaction energy U increases algebraically 
with the lattice depth. There is another important pa- 
rameter, the number of particles per site v (also called as 
filling factor) that is implicit in Eq. ([T]). The ratio U/ {vJ) 
controls the ground state phase of the BH model. When 
U /{vJ) < 1, the SF phase is favored at any filling factors. 
When U /{vJ) increases at an integer filling factor such 
that U /{vJ) > 1, quantum fluctuations drive the system 
to the MI phase. 

Since the SF to MI transition is one of the most re- 
markable phenomena emerging in the BH model, many 
previous studies have made efforts to determine the tran- 
sition point both numerically [l^ - [20| and experimen- 
tally 0, H, H Eli • Especially, theorists have accurately 
calculated the ratio ([// J)c at the transition point for dif- 
ferent filling factors and dimensionalities. In 2D and 3D, 



quantum Monte Carlo simulations have provided (U/ J)c 
at v = 1 [3, [l^ while the transition points at arlDitrary 
filling factors have been calculated with use of the strong- 
couphng expansion (SCE) techniques [H, [13 • In ID, 
{U / J)c has been determined at low filling factors, namely 
at = 1 [ll-i^l and = 2 di, iOl, using the quasi- 
exact numerical methods of density-matrix renormaliza- 
tion group (DMRG) [2l| and time-evolving block decima- 
tion (TEBD) [l^l. However, the transition points in the 
high filling region [v > 3) are yet to be obtained because 
of higher computational cost. Difficulty in treating the 
high filling region stems also from the fact that SCE fails 
to give an accurate transition point in ID because the 
transition is of the Berezinski-Kosterlitz-Thouless (BKT) 
type. 



In this paper, we calculate the critical points of the SF 
to MI transition in ID for arbitrary integer filling factors, 
including the region of very high filling factors in which 
the BH model is equivalent to the 0{2) quantum rotor 
model [1^. We emphasize that the determination of the 
transition points in the quantum rotor regime is impor- 
tant in the sense that the quantum rotor model effectively 
describes a regular array of Josephson junctions [23 | and 
liquid '''He absorbed in nanopores [25l - [27| . In order to 
locate the transition point, one usually calculates static 
quantities such as the single-particle density matrix (b\bi) 
and the density-density correlation function {fijfii) |18l - 
[20| . In contrast, here we use a dynamic quantity, that 
is, energy splitting A caused by tunneling between two 
states with macroscopically distinct currents (or winding 
numbers) p8| . In our previous work [2^ , we have shown 
that the region of high filling factors can be efficiently 
treated with TEBD by imposing a lower bound for the 
occupation number at each site in addition to an upper 
bound. It has been also shown that the tunneling split- 
ting can be accurately computed by means of TEBD. 
Using these prescriptions, we calculate A to determine 
the transition points. We find that the transition point 
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function of v is well approximated by 



described as 
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where D denotes the dimensionality of the system (e.g. 
D ~ 1 for ID), and the constants a, 6, and c arc nu- 
merically determined. We show that one can express 
(U/ (Di/J))r. also in 2D and 3D, which has been obtained 
in Ref. [ij, in the form of Eq. ([2]) with different values of 
the constants. Given these results, one can immediately 
know the transition points for any filling factors and any 
realistic dimensions. 

This paper is organized as follows. In Sec. HIl we in- 
troduce the system and the model considered here and 
explain how to determine the transition point from the 
tunneling splitting. Wc qualitatively discuss the relation 
between the tunneling splitting and the SF to MI tran- 
sition from a perspective of the instanton method. In 
Sec. IIII[ in order to examine the performance of the sug- 
gested procedure, we apply it to the ID hardcore BH 
model with nearest-neighbor interactions that is exactly 
solvable by means of the Bethe ansatz. In Sec. lIVi we cal- 
culate the transition point in the BH model of Eq. ([ij as 
a function of the filling factor. In Sec. |Vl we summarize 
the results. 



Here v is the sound speed. 

It is well-known that the SF phase is favored when 
K > 2/nl and that the SF to MI transition of the BKT 
type occurs at K = 2/n^ [s^. At the transition point, 
the renormalization of the TL parameter leads to a log- 
arithmic correction on the scaling of A as [2^ 
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In the next section, we will see that the inclusion of the 
logarithmic correction is important to accurately calcu- 
late the transition point. 

Given the scaling formulas of Eqs. ([3]) and ([5]), one can 
determine the transition point from the energy splitting 
as follows. First, using TEBD, we numerically compute 
A as a function of L in the way described in Ref. [l^ We 
next fit a function 



-B 



fiL) 



AL 
InL 



(6) 



II. HOW TO DETERMINE THE TRANSITION 
POINT FROM THE ENERGY SPLITTING 

Wc consider a system of ID lattice bosons in a ring- 
shaped geometry, i.e. with a periodic boundary. We as- 
sume for the moment that the system is in the SF phase. 
In the classical limit, which corresponds to the limit of 
U/{vJ) — >■ in the BH model, a state with a finite ho- 
mogeneous current can be metastable, and its quasimo- 
mcntum per particle is discretizcd as p = 2n7rft/(L(i), 
where the integer n is the winding number, L the num- 
ber of lattice sites, and d the lattice spacing. We sup- 
pose a situation that two states with different winding 
numbers, say ni and n2, are degenerate. We define the 
winding number difference between the two states as 
rid = \ni — n2\- When n^iV is an integer and U/{i>J) 
is finite, Umklapp scattering processes can induce phase 
slips via quantum tunneling to couple these two states. 
As a result, the degeneracy is broken and there emerges 
the energy splitting A that quantifies the tunneling rate. 
Applying the instanton techniques [13, [Hj to the phe- 
nomenological Tomonaga-Luttinger (TL) liquid model, 
Kashurnikov et al. have derived the following scaling for- 
mula of the energy splitting at the smallest possible wind- 
ing number difference rid with respect to L psl ]. 

A (X L-"d^+i (3) 

where the TL parameter K is defined such that the ef- 
fective action for the phase of the bosonic field 9{x, r) is 



to the numerical data, where A and B are free param- 
eters, and extract the exponent B. We determine the 
transition point from the condition that B = 1. 

Let us explain a reason why the exponent B is equal to 
unity at the transition point from a different viewpoint, 
which is the relation between the SF to MI transition 
and the validity of the instanton formula of Ec^. ([3|) . An 
important point is that the instanton techniques used 
to derive Eq. Q are based on the so-called dilute gas 
approximation (DGA), in which the instantons are as- 
sumed to be well separated from each other in the path- 
integral trajectories that contribute to the partition func- 
tion [13, [m. In our system of ID lattice bosons, since 
an instanton is regarded as a vortex in the space-time 
plane, the breakdown of DGA means that many vortices 
are created in the space-time coordinate so that they de- 
stroy the long-range coherence of bosonic phases, leading 
to the quantum phase transition to the Mott insulator. 
In short, the breakdown of DGA signals the Mott tran- 
sition [33|. In general, the condition under which the 
dilute gas approximation is valid is that the size of an in- 
stanton along the imaginary-time axis ti is much smaller 
than the tunneling time h/ A. In the present case, since 
the winding number difference is of the order of unity, 
the instanton size is of the order of the system size, i.e. 
Ti Qc L. Meanwhile, A is given by Eq. Obviously, 
when B > 1, the condition that ri ^ ft/A is held in the 
thermodynamics limit such that DGA is valid, i.e. the 
system is in the SF phase. Thus, the condition for the 
Mott transition is given by _B = 1. 
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FIG. 1: (color online) The time evolution of the current ve- 
locity v{t) in the dynamics of the model of Eq. 0, where 
L = 40 and V/J = 1.8. 
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FIG. 2: (color online) The time evolution of the overlaps Oi (t) 
(blue solid line) and 0-i{t) (black dashed line) in the dynam- 
ics of the model of Eq. 0, where L = 40 and V/J= 1.8. 



III. HARDCORE BOSE-HUBBARD MODEL 
WITH THE NEAREST-NEIGHBOR 
INTERACTIONS 

In this section, we analyze the ID hardeore BH model 
with the nearest-neighbor interactions, 

L L 

in order to illustrate that the critical point of the ID 
superfluid to insulator transition can be accurately cal- 
culated along the procedure described in the previous 
section. Here, V is the nearest-neighbor interaction, and 
Cj and TTij are the creation and number operators of a 
hardcore boson at site j. We also include the phase twist 
e~'^ in the hopping term in order to control the winding 
number of states. The hardcore constraint means that 
the maximum occupation number at each site is unity. 
This constraint leads to the identities between the oper- 
ators of hardcore bosons and those of 1/2-spins, namely 
Sj = rhj — 1/2 and S~ = cj, which tell us that the model 
of Eq. ([7]) is equivalent to the spin- 1/2 XXZ model and is 
exactly solvable by means of the Bethe ansatz [s^]. Ac- 
cording to the exact solution, there is a quantum phase 
transition between the SF and the density- wave insulator 
at half filling due to the competition between V and J. 
When V/J = 0, the particles favor to be delocalized and 
the system is in the SF phase. When V/J increases, the 
Umklapp scattering tends to localize the particles, and 
the transition to the insulating state occurs at V/ J = 2. 
In the following, we show that our procedure provides a 
numerical value of the transition point that is very close 
to the exact one. 

Since the SF to MI transition occurs at = 1/2 in the 
model of Eq. the minimum winding number differ- 
ence in possible phase slip processes is rid = 2. In order 
to obtain the energy splitting for = 2, we first pre- 
pare the ground state of Eq. ([7]) with 9 = 2tt/L, which 
is a state with winding number n = 1. For dealing with 
our system with a ring-shape geometry, we use TEBD 
for a periodic boundary condition [35| . Taking the state 



with n = 1 as the initial state and setting 6* = 0, we 
next carry out the real-time evolution. Since two states 
with n = 1 and n = — 1 are degenerate when 6 = 0, 
it is expected that the supercurrent dynamics exhibit a 
coherent oscillation between these two states induced by 
quantum tunneling. To demonstrate this, we calculate 
the time evolution of the current velocity given by 

^' = ^E<4''j+i-'^-c-)' (8) 

j 

where N is the total number of particles. In Fig. [I] v{t) 
for L = 40 and V/J = 1.8 is shown. We clearly see that 
the velocity oscillates between v{t = 0) and —v{t = 0), 
i.e. between the states with n = 1 and n = —1. We 
also calculate the overlap 0„(t) = |('I'„|^'(^))|^ of the 
wave function with the ground state of the Hamil- 
tonian ([7]) with 6 = 2'im/ L. In Fig. [2l we show overlaps 
Oi{t) and O- i{t), which reconfirm that the wave function 
4'(i)) coherently oscillates between |$i) and l^-i). We 
note that similar quantum-tunneling dynamics have been 
found also for quantum vortices in anisotropic traps f36j 
and supercurrents in two-color optical lattices [37| . 

Once the tunneling dynamics in real time are obtained, 
we can extract the energy splitting A by fitting Oi{t) 
using the function 

.9(t)-Ccos2 (^^i^ +F (9) 

where A, C, and F are free parameters. 

In Fig.[3l we plot the energy splitting between |<I>i) and 
as a function of the number of lattice sites L up 
to L = 100. The numerical data are very well fitted by 
the function of Eq. ([6]). From this fitting, we extract the 
exponent B and calculate it varying V/ J, as shown by 
the red squares in Fig. As expected, B monotonically 
decreases with V/ J . The condition that B = 1 gives the 
transition point as {V/J)c = 1.952 ± 0.008. This value 
is so close to the exact result, (V^/J)c, exact = 2, that 
the relative error, \ {V/J)c - (X^/J)c,oxact|/(V"/ J)c, exact, is 
smaller than 3%. 

We emphasize that our procedure provides an accurate 
value of the transition point even when the system size is 
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FIG. 3: (color online) The energy splitting versus the number 
of lattice sites L for the ID hardcore Bose-Hubbard model 
of Eq. (0 at half filling. We take VjJ = 1.8 (black dia- 
monds), 2.0 (blue squares), and 2.2 (red diamonds). The 
solid lines represent the best fits to the respective numerical 
date with the function of Eq. (O, where {A, B) = (2.80, 1.11), 
(2.28,0.965), and (1.84,0.823) for V/J = 1.8, 2.0, and 2.2. 



TABLE I: The critical values {U/{uJ))c of the ID BH model 
for several values of the filling factor p. 




FIG. 4: (color online) The exponent B versus V/J for the 
ID hardcore Bose-Hubbard model of Eq. ((Tjl at half filling. 
B is extracted by fitting the numerical data of A versus L 
upto L = 100. We use the fitting function of Eq. ((6} with the 
logarithmic correction for red squares, while we use Eq. (|10p 
without the correction for the black circles. The solid lines 
are guides to the eyes. The dashed line represents B — 1, 
which is the condition to determine the transition point. 
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FIG. 5: (color online) The red circles represent the critical 
point (U/(yJ))c versus the filling factor. The solid line rep- 
resents the best fit to the numerical data with the function 
of Eq. ((2)1, where the numerical constants are obtained as 
(a, 6, c) = (2.16, 0.97, 2.13). The transverse axis is depicted in 
a logarithmic scale. Notice that the critical point exists only 
at integer fillings although the solid line is continuous. In the 
inset, the same data are presented together with the transi- 
tion point calculated with the analytical formula of Eq. (|12p 
(green dashed line). 



relatively small. For instance, when we take the number 
of lattice sites up to i = 48 for the fitting of the energy 
splitting, we find that iy/J)c = 1.925 ± 0.013 and that 
the relative error is still within 5%. This is a clear advan- 
tage of the use of the energy splitting over the correlation 
functions that require a substantially larger system size. 
We also note a disadvantage of our procedure that one 
needs a system with periodic boundaries in which the bi- 
partite entanglement entropy is twice as large as that in 
a system with open boundaries. 

In order to corroborate the necessity of the logarithmic 
correction in the fitting function of Eq. we instead 
use another function that docs not include the correction, 

J{L) = AL-^, (10) 

for the fitting. The black circles in Fig. 2] represent B 
extracted with f{L), and the value of the transition point 



in this case is (V/J)c ~ 2.37 ± 0.02. The relative error 
is ^ 20%, which is much larger than the case with the 
logarithmic correction. 



IV. THE BOSE-HUBBARD MODEL AT 
ARBITRARY INTEGER FILLINGS 

In this section, we determine the critical point of the 
SF to MI transition of the ID BH model (P for arbitrary 
integer fillings, which is the main goal of the present pa- 
per. For this purpose, we start with the ID BH model 
with a phase twist, 

L L 
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FIG. 6: (color online) The critical value {U / {DvJ))c versus 
the filling factor v for 2D (blue squares) and 3D (black dia- 
monds). The solid lines represent the best fit to the numerical 
data with the function of Eq. The dashed lines represent 
the transition points calculated from the analytical formula 
of Eq. (|12p . The transverse axis is depicted in a logarith- 
mic scale. The numerical data for 2D and 3D are taken from 
Ref. ITtI . Notice that the critical point exists only at integer 
fillings although the solid line is continuous. 



In the case of integer fillings, the minimum winding num- 
ber difference in possible phase-slip processes is = f . 
In our previous work, we have shown a way to obtain the 
energy splitting for rid = 1 [291 using TEBD, which we 
follow in the present paper as well. We first calculate the 
ground state of Eq. ([Tl]) with = 2'k/L to obtain a state 
with n = 1. Taking this state as an initial state and set- 
ting 9 = tt/L, we calculate the real-time evolution. Since 
|$i) is degenerate with |<i>o) at 9 = n/L, the dynamics 
exhibit a coherent oscillation between these two states. 
We extract the energy splitting from the oscillation, and 
calculate it as a function of the number of lattice sites up 
to L = 48. The transition point is determined from the 
L-dcpendcnce of A as done in the previous section. 

In Table U and Fig. [5] (red circles), we show the critical 
point {U/{vJ))c as a function of v. Large scale DMRG 
analyses by Ejima et al. have provided the the most re- 
cent benchmark values of the critical points ioi v = 1 
and 2 as {U/{vJ))c = 3.28 and 2.78 Our results 

agree well with them to the extent that the relative devi- 
ation is within 5%. Since the ratio U / [vJ) quantifies the 
strength of quantum fluctuations in the quantum rotor 
limit {v oo) [1^, {U/{vJ))c is expected to be con- 
verged to a certain value when v is sufficiently large. In- 
deed, the critical value monotonically decreases with v 
and approaches an asymptotic value when v ^ 1. At 
V = 1000, {U/{vJ))c is well converged to the asymptotic 
value that corresponds to the quantum rotor limit. 

In Ref. [TtI . it has been shown that an analytical formula 
accurately approximates the transition point for D > 2 



Dimensionality: D 


(a,b,c) 


1 
2 
3 


(2.16, 0.97, 2.13) 
(5.80, 2.66, 2.19) 
(6.70, 3.08, 2.18) 



TABLE II: The constants in the function of Eq. ((2]) obtained 
by the best fit to the numerical data for ID, 2D, and 3D. 
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(13) 



is the transition point obtained by a mean-field the- 
ory [l3]. Although it has not been argued that this for- 
mula is valid for ID, it is worth checking whether it works 
in ID or not. The green dashed line in the inset of Fig. [5] 
represents the transition point given by Eq. (|12p . Obvi- 
ously, the formula fails; the relative error is almost 100%. 
This is not totally unexpected because the transition in 
ID is special in the sense that it is of the BKT type while 
the transitions in higher dimensions are of the second- 
order. Instead of Eq. ([T^ . we show that another ana- 
lytical formula of Eq. ([2]) well approximates {U/{vJ))c 
versus v. As seen in Fig. [Sj the fitting with the function 
of Eq. ([2]) agrees with the data to the extent that the 
deviations are within the size of the data points. 

The function of Eq. ([2]) is a good approximation also 
for the transition points in 2D and 3D. To show this, we 
depict in Fig.[6]the numerical data of [U /{DvJ))c versus 
V for 2D (blue squares) and 3D (black diamonds), which 
are obtained using SCE in Ref. [13, together with the best 
fit with the function of Eq. ^ represented by the solid 
lines. The transition points calculated from Eq. ([T2|) are 
also plotted with the dashed lines. There we see that the 
function of Eq. ^ approximates the transition points 
as accurately as Eq. (IT^ . The values of the numerical 
constants (a, 6, c) in the fitting function are summarized 
in Table HH 



V. CONCLUSIONS 

In summary, we have studied the superfiuid to Mott 
insulator transition in a system of one-dimensional lat- 
tice bosons. We have shown that the transition points in 
ID can be accurately determined from the energy split- 
ting between two degenerate states with distinct wind- 
ing numbers. We obtained the transition points for the 
Bose-Hubbard model with arbitrary integer filling fac- 
tors, including the high filling limit corresponding to the 
quantum rotor regime. We have found a simple analyt- 
ical formula of Eq. ^ that well approximates the tran- 
sition point versus the filling factor for ID, 2D, and 3D. 
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With this formula, the transition points can be obtained 
easily and immediately for any fillings and any realistic 
dimensionalities . 
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